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ABSTRACT 

We critically examine how the evolution of the matter density field in cosmolog- 
ical simulations is affected by details of setting up initial conditions. We show that 
it is non-trivial to realise an initial distribution of matter in TV-body /hydrodynamic 
simulations so that the baryon and dark matter density fluctuations and their veloc- 
ities evolve consistently as theoretically predicted. We perform a set of cosmological 
simulations and use them to distinguish and verify an appropriate method for generat- 
ing initial conditions. We show that a straightforward way of applying the Zel'dovich 
approximation to each component using distinct transfer functions results in an incor- 
rect growth of density fluctuations and that it is necessary to correct velocities at the 
initial epoch. The unperturbed uniform particle distribution must be also generated 
appropriately to avoid tight coupling of the baryonic and dark matter components. 
We recommend using independent "glass" particle distributions, using distinct trans- 
fer functions for baryons and dark matter, and taking into account the difference in 
the velocity fields at the initialisation epoch. The proposed method will be useful for 
studies of the evolution of the intergalactic medium and the formation of the first 
cosmological objects using numerical simulations. 
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1 INTRODUCTION 

High-precision measurements of the Cosmic Microwave 
Background (CMB) radiation by the WMAP satellite pro- 
vide strong support for the so-called standard model of cos- 
mic structure formation, in which Cold Dark Matter (CDM) 
and dark energy dominate the Universe. The matter distri- 
bution at the decoupling epoch has been directly probed by 
measurements of CMB anisotropies, and it has been con- 
firmed that the observed matter distribution is consistent 
with the predictions of popular inflationary theories (Spergel 
et al. 2003). Within this theoretical framework, structure 
formation in the early Universe is described by gravita- 
tional amplification of small initial density fluctuations; non- 
linear evolution leads to the formation of dark matter ha- 
los, followed by hydrodynamic processes such as gas infall 
into gravitational potential wells, shock heating, and radia- 
tive cooling (e.g. Yoshida et al. 2003). Precise modelling of 
the angular power spectrum of the CMB anisotropies deter- 
mines the relative densities of the baryonic and non-baryonic 
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components. Furthermore, accurate solutions of the non- 
linear Boltzmann equations (Hu & Sugiyama 1995; Ma & 
Bertschinger 1995; Seljak & Zaldarriaga 1996) show that 
there is a substantial difference in the distribution of baryons 
and dark matter at the decoupling epoch. Indeed, recent 
detailed numerical studies (Yamamoto, Sugiyama & Sato 
1998; Singh & Ma 2002) show that the differences remain 
until much lower redshift z ~ 20 for high baryon densities. 
In fact, the ratio of baryon density to CDM density deter- 
mined by the WMAP data, (Spergel et al. 2003) is as large as 
fibaryon/ficdm ~ ^-2. Thus these results clearly contradicts 
the assumption, often made in the context of structure for- 
mation, that baryons trace the dark matter; i.e. that local 
baryon densities are just proportional to those of the dark 
matter, and their velocities are identical both in direction 
and in amplitude. 

Numerical simulation of structure formation is aimed at 
revealing how the large-scale features of the Universe formed 
and evolved from high to low redshifts. Conventionally, cos- 
mological iV-body simulations are started from an initial 
matter distribution - a random Gaussian realisation - gen- 
erated using an input power spectrum for the total matter 
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density. Although this is obviously the most plausible way 
of generating initial conditions for iV-body simulations with 
only a single collisionless component, it is non-trivial to re- 
alise the matter distribution in simulations that employ two 
(or more) components so that the density and velocity fields 
of both components evolve consistently. A common belief is 
that details in the initial configuration are erased by non- 
linear evolution and hence do not matter if one is interested 
in the structure of non- linear objects (i.e. "dark halos") at 
low redshift. However, this reasoning may not apply at high 
redshifts or for the intergalactic medium (IGM) where the 
density fluctuations on relevant scales are in the linear to 
mildly non-linear regime. 

In numerical studies of the properties of the IGM, such 
as those concerned with the reionisation history of the Uni- 
verse (e.g. Gnedin & Ostriker 1997; Gnedin 2000; Sokasian 
et al. 2002, 2003), it is important to accurately compute the 
thermal properties of baryons in low density regions. Various 
analytic methods have been developed to obtain the local 
densities and temperatures of the IGM from the distribution 
of dark matter that can be computed relatively accurately 
by either direct iV-body simulations or linear theory (e.g. Bi, 
Borner & Chu 1992; Jones 1996; Hui & Gnedin 1997; Nusser 
2000; Matarrese & Mohayaee 2002). On the other hand, in 
cosmological iV-body/hydrodynamic simulations, simplify- 
ing assumptions are often made. For example, in cosmolog- 
ical Smoothed Particle Hydrodynamics (SPH) simulations, 
gas particles are either put on top of dark matter particles, 
or simply displaced by half the mean interparticle separation 
before perturbations are imposed (e.g. Katz et al. 1996). A 
naive way of improving upon this simple approach would 
be to apply the Zel'dovich approximation to baryonic and 
dark matter components separately using separate transfer 
functions for each. It is unclear, however, if the generated 
initial particle distributions and the velocity field reproduce 
consistent results with predictions of full non-linear theory. 

Another closely related issue is that, in particle simula- 
tions, a smooth density field must be represented by discrete 
mass elements. It is well-known that discretisation itself im- 
poses some limitations on the accuracy of numerical simula- 
tions (Splinter et al. 1998; Hamana, Yoshida & Suto 2001; 
Baertschiger et al. 2002). Baertschiger & Sylos-Labini (2001) 
show that the discreteness of the particles contributes to the 
evolution of power-law clustering on small scales. Gotz & 
Sommer-Larsen (2002) argue that even small differences in 
unperturbed particle distributions can yield quite different, 
in some cases unfavourable, results in the nonlinear evolu- 
tion of systems (see their Fig. 1). 

These facts and questions motivate the need for a sys- 
tematic study in order to develop and verify a proper method 
for setting up initial conditions for multi-component cosmo- 
logical simulations. In the present paper, we examine how 
details in the initial set-up affect the evolution of matter den- 
sity fluctuations. We specifically study the power spectra of 
baryonic and dark matter components and the abundance 
of dark matter halos using a set of iV-body simulations. We 
show that differences in initial configurations indeed result 
in substantial differences in the evolved density field and its 
growth rate. We compare the numerical results with theo- 
retical predictions and determine an improved method for 
generating initial conditions. 

The rest of the paper is organised as follows. In Section 
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Figure 1. The density and velocity power spectra for CDM (solid 
line) and for baryons (dashed line) at z = 100 for a standard 
ACDM model with C cdm = 0.26, O baryon = 0.04, fl A = 0.7, H = 
70 km s — x Mpc — x . In the top panel, the relative amplitudes of 
the density to velocity power spectra are scaled arbitrarily. The 
middle panel shows the ratio of the two density power spectra and 
the bottom panel shows that of the two velocity power spectra. 



2, we present the basic theory of the evolution of density 
fluctuations in the early Universe. In Section 3, we describe 
the simulation set. We show our numerical results in section 
4. Discussion is given in Section 5. 



© 0000 RAS, MNRAS 000, 000-000 



Baryon density fluctuations in cosmological simulations 3 



2 THEORY 

2.1 Evolution of baryon density fluctuations 

The evolution of the cosmological baryon perturbations be- 
fore the decoupling epoch is described by a number of physi- 
cal processes which couple the perturbation fields of baryons, 
CDM, photons and neutrinos in a complex manner. One thus 
needs to solve the full set of coupled nonlinear Einstein- 
Boltzmann equations in order to obtain an accurate re- 
sult. Recent detailed studies by Liu et al. (2001) and Singh 
& Ma (2002) show that, although density fluctuations on 
scales smaller than the photon diffusion scale (Silk 1968) are 
strongly suppressed, the growth of the baryon density fluc- 
tuations are eventually accelerated during the drag epoch 
because the photon-baryon coupling breaks down by recom- 
bination. It is also important to note that second-order ef- 
fects arising from the coupling of the velocity fields to the 
baryon density field can enhance the growth of small scale 
fluctuations (Shaviv 1998; Liu et al. 2001; Singh & Ma 2002), 
although this effect plays an important role only on very 
small length scales (k > 1000Mpc _1 ). 

After recombination, gravitational infall causes the 
baryon density fluctuations to catch up with those of the 
dark matter. Yamamoto et al. (1998) show, however, that 
the growth of perturbation modes on length scales smaller 
than the Jeans length is delayed due to the initial oscillatory 
behaviour of the baryons, and thus the residual differences in 
the fluctuation amplitudes between baryons and dark mat- 
ter remain until redshifts z ~ 20— 100. This is clearly shown 
in Figure 1. There we plot the density power spectra and the 
velocity power spectra for baryons and CDM for a conven- 
tional ACDM universe (parameters as summarised in Sec- 
tion 3) at z = 100. The power spectra are computed us- 
ing the Boltzmann solver of Sugiyama (1995), as described 
in Yamamoto, Sugiyama & Sato (1998). The density power 
spectra for baryons and CDM differ by about a factor of two 
for 1 < k < 10, whereas the difference in the velocity power 
spectra is less than 10 percent on these scales (compare the 
middle panel with the bottom panel) . 



2.2 The Zel'dovich approximation 

The standard procedure for setting up initial conditions 
for cosmological JV-body simulations is to perturb a uni- 
form particle distribution using the Zel'dovich approxima- 
tion (Zel'dovich 1970). An extensive study of the applicabil- 
ity of the Zel'dovich approximation is presented in a series 
of papers by Valageas (2002a, b and references therein). We 
summarise the elements of the method in this section. 

Let the initial Lagrangian coordinate of a particle in an 
unperturbed distribution be q. Then each particle is subject 
to a displacement corresponding to a density perturbation. 
In the Zel'dovich approximation the Eulerian coordinate of 
the particle at time t is 



Table 1. Simulation sets 



r(t,q)=a(t)[q-6(t)V q *o(q)], 



(1) 



where r = a(t)x with x being a comoving coordinate and 
b(t) is the growth factor. The velocity field is obtained from 
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1 Transfer function for the total matter 

The quantity $o(q) is related to the density perturbation S 
in the linear regime via the Poisson equation: 



bVq$o. 



(3) 



v = = — o6V q $o(q)- 



(2) 



It is important to point out that the initial velocity field and 
the displacement field are linearly related via equation (2). 
Namely, particles in the Zel'dovich approximation execute 
motion on straight line trajectories. Therefore, if one uses 
different transfer functions for two components (e.g. CDM 
and baryons) while keeping the phase information for a ran- 
dom Gaussian field identical, the generated displacement 
fields and hence the velocity fields for the two components 
should differ by an amount corresponding to the difference 
in the amplitudes of the transfer functions. However, after 
the decoupling epoch, baryons fall into gravitational poten- 
tial wells that are dominated by the dark matter, whereas 
the growth of dark matter density fluctuations is slightly 
delayed due to the baryonic matter which is more homo- 
geneously distributed. Then, the difference in velocities of 
baryons and dark matter decays quickly after the decoupling 
epoch (compare the middle panel and the bottom panel in 
Figure 1). 

Nusser (2000) and Matarrese & Mohayaee (2002) re- 
cently derived analytic solutions for the evolution of linear 
density density perturbations for both baryons and dark 
matter. While their solutions fully describe the evolution 
of linear density fluctuations of a two-component fluid for a 
few particular cases, they are not directly applicable to gen- 
crating initial conditions for direct iV-body simulations in 
more general cases. In the present paper we develop a more 
practical method. The key question is how can we generate 
initial conditions using solutions of a full non-linear calcu- 
lation ? 

We test five different methods for setting up initial con- 
ditions in order to distinguish a correct method (if one ex- 
ists) and to quantify its accuracy. We describes the technical 
details in the next section. 



3 THE A-BODY/SPH SIMULATIONS 

We use the parallel iV-body/SPH code GADGET (Springel, 
Yoshida & White 2001), in its "conservative entropy" ver- 
sion (Springel & Hernquist 2002, 2003). All our simulations 
employ 2 x 128 3 particles (CDM and non-radiative gas com- 
ponents) in a cosmological box of 4 fe _1 Mpc on a side. We 
work with a flat A-dominated Cold Dark Matter universe, 
with matter density f2 cdm = 0.26, Qb = 0.04, cosmological 
constant Qa = 0.7 and expansion rate at the present epoch 



© 0000 RAS, MNRAS 000, 000-000 



4 N. Yoshida, N. Sugiyama & L. Hernquist 



Ho — 70km s _1 Mpc _1 . The qualitative features of our re- 
sults are not affected by the adopted cosmology, although 
the degree of the numerical artifacts we discuss depend on 
the ratio of f2 c d m /^b- Table 1 summarises the five methods 
we consider for the set-up of the initial conditions. We de- 
scribe each specifically in the following subsections. We note 
that for all the simulations we set the parameters governing 
the accuracy of the force calculation and time stepping to be 
very conservative, following Springel et al. (2001) and Power 
et al. (2003). We set the gravitational softening length e ao f t 
= 3 h~ kpc ~ one tenth of the mean inter-particle separa- 
tion. 



3.1 Unperturbed particle distribution 

There are two popular methods for distributing particles 
"uniformly". One is to place particles on a rectangular 
grid. While used most often, a grid distribution has the 
unfavourable feature of imposing a particular direction in 
the particle distribution. In other words, a grid distribu- 
tion is not isotropic, although it is homogeneous on scales 
larger than the mean inter-particle separation (see e.g. 
Baertschiger & Sylos-Labini 2001). G6tz & Sommer-Larsen 
(2002) report that grid distributions cause peculiar problems 
when used for warm dark matter simulations such as those 
of Bode, Ostriker & Turok (2001). Although the degree of 
these problems is somewhat uncertain and is likely to be case 
dependent, it is preferable to avoid these issues altogether by 
using "glass" particle distributions. Following White (1994), 
we generate glass initial conditions by evolving a gravita- 
tionally interacting iV-particle system in an expanding back- 
ground with the sign of gravity reversed such that each par- 
ticle feels repulsive forces from all the other particles. When 
the system reaches a quasi-equilibrium state, unphysical 
clumping in the initial distribution is substantially damped, 
while maintaining uniformity. The particle distribution is 
also nearly isotropic and hence does not possess any pre- 
ferred direction. 

For all the numerical experiments in the present paper, 
we use glass particle distributions. We generate two inde- 
pendent glass distributions and use them for gas and CDM 
particles in Runs B, C, and E, whereas we use identical glass 
particle distributions for both components in Run A and D 
(see Table 1). 

3.2 Transfer functions 

In the context of structure formation in CDM models, the 
density power spectrum at a given epoch is related to the 
transfer function by 

P{k) = Ak n T 2 (k), (4) 

where A is the normalization factor. We consider the 
Harrison-Zel'dovich primordial power spectrum with n = 1. 
Figure 1 shows the power spectrum 4ivk 3 P(k) at z- ln i t = 100. 
On very large length scales (k < 0.01/i/Mpc), the power 
spectra for baryons and CDM have the same amplitudes, 
whereas on intermediate length scales (1 < k < 100) the 
amplitudes for baryon density perturbations are smaller by 
about a factor of two at z = 100. Also the power spectrum 
for baryons shows characteristic wiggles at 0.01 < k < 0.2 



which is induced by Jeans oscillations before the recombi- 
nation epoch. 

For Runs A, B, and C, transfer functions are given sep- 
arately for the CDM and the baryonic components, whereas 
for Run D, we use a single transfer function computed for 
the total matter density fluctuations. The matter transfer 
function for Run D is computed from 

T(k) = ^T b (k) + ^T cdm (fc). (5) 

"total "total 

3.3 Phases 

In standard inflationary theories that predict adiabatic-type 
density fluctuations, the distribution of various forms of 
matter is described by essentially a single scalar function, 
and thus the phases for baryons and CDM are identical for 
all the perturbation modes. Although we have no compelling 
reason to assume the phases are different on any length scale, 
it is interesting to examine the effect of varying the phases. 
In order to make comparisons, we assign random, indepen- 
dent phase information for baryons and CDM in Run E. 
Other configurations for Run E are identical to those of Run 
B. We compare these simulations in Section 5. 

3.4 Velocity power spectrum 

As we have discussed in Section 2, the differences in the ve- 
locity power spectra of baryons and CDM are smaller than 
the differences in fractional amplitudes of density fluctua- 
tions (Figure 1). Since the fractional difference is within 10% 
and nearly a constant at Zi n i t = 100 over most of the relevant 
scales, we approximate the velocity fields of both CDM and 
baryons by the velocity field of the total matter. In practice, 
we compute the velocity field using the transfer function for 
the total matter while keeping the phase information iden- 
tical to that used in generating the spatial displacements. 
Note that this is equivalent to multiplying the velocity term 
in equation (2) by a scalar factor. The scalar factor could, 
in principle, be made scale-dependent, as is done by Klypin 
et al. (1997) for Mixed Dark Matter simulations. 

4 RESULTS 

Figure 2 shows the power spectra of the initial mass distri- 
bution for the four simulations (Run A, B, C, D), and the 
subsequent evolution is shown in Figure 2 (for z = 50) and 
Figure 3 (z = 20, 30). The filled circles are measured power 
spectrum for CDM and the open diamonds are for baryons. 
Note that the particular random realization for these runs 
shows a dip at k ~ 3ft/Mpc. This is simply due to ran- 
domization when assigning the fluctuation amplitudes for 
each perturbation mode, and irrelevant to the analysis in 
the present paper. Theoretical predictions are shown by the 
solid line (CDM) and by the dashed line (baryon) for each 
output redshift. The measured power spectra at z = 100 
show a steep turn-over at large wavenumbers k > 20 due to 
the Poisson noise which scales with the particle number N 
as (\5 k \ 2 ) = 1/N. (Note that we plot 47rfc 3 P(fc) in Figures 2 
and 3.) Although the level of the Poisson noise could be re- 
duced by using a larger number of particles, we refrain from 
carrying out such large, costly simulations. Since we focus 



© 0000 RAS, MNRAS 000, 000-000 



Baryon density fluctuations in cosmological simulations 5 



0.1000 



0.0100 



if. 



0.0010 



0.0001 



■ 1.4 
ffi 1.2 



0.1000 



0.0100 



0.0010 



0.0001 



: Run A 








""^« - ' ' 


'"^ 




' — • 

,-'0 






-''o 


CDM : 


.,-''z'=100 

iii 




baryon 


1 


k [h/Mpc] 


10 


r z=50 ° 


o o ° 
• »....« 


°0<XK> ^ 


• 




.^y 


1 


k [h/Mpc] 


10 


- i i i i i 

; Run C 




11 i 




, - O' 

--'© 




\, - ' ' 


,-'0 


^<y - - - 




-''o 


CDM 


,-'z=100 




baryon 


1 


k [h/Mpc] 


10 



0.1000 



0.0100 



Q- 



0.0010 



0.0001 



-rl 1 - 4H 

S 1.2H 

Q. 1.0 L 

§ 0.8 r 

=5 0.6^ 



0.1000 




k [h/Mpc] 



z=50 



• • ••••••• 



1 10 
k [h/Mpc] 



0.0100 



CL 

Si 



0.0010 



0.0001 



- J 

; Run D 












: z=50^<^ 


• ^^-^ 








-"' ♦ 








CDM : 


- z=100 


baryon ■ 




total 



10 



k [h/Mpc] 



"D 

=1 



1.4 r 
1.2 H 
1.0 r- H 
0.8 r 
0.6 H 



z=50 



=5 



1.4 r 

1.2 h 
1 .0 

0.8 r 
0.6 r 



10 



k [h/Mpc] 



z=50 



♦ ♦ ♦ 



10 



k [h/Mpc] 



Figure 2. The measured power spectra for the baryonic and CDM components are compared with the theoretical predictions at the 
initial epoch (z = 100) and at z = 50. In each panel the top portion shows the power spectra 4irk 3 P(k) and the bottom portion shows 
the measured growth rate normalised by the theoretically expected value. Note the vertical axis in the bottom portion is in linear scale. 
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Figure 3. As for Figure 2, but for z = 20 and z = 30. The power spectra at z = 20 are scaled vertically by a constant factor of five, in 
order to clarify the plot 
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on density fluctuations on large length scales, k < 20, where 
the power spectra and the growth rates can be measured ro- 
bustly, the shot noise effect is not important in the analysis 
presented in the following. 

In Run D, baryons and CDM have exactly the same ve- 
locities both in amplitude and direction. For this simulation 
the gas particles are simply put on top of the CDM parti- 
cles and thus they move together initially. Since we used the 
transfer function for the total matter density in Run D, the 
measured power spectrum for CDM is smaller than the the- 
oretical prediction (solid line), whereas that for the baryons 
is larger than the theoretical prediction (dashed line). These 
differences can be easily accounted for by the ratio of fibaryon 
and f2 c dm to fltotai' 

In the bottom portion of each panel in Figure 2, 3, we 
plot the growth rate of the densities. We measure the growth 
rate between Zmit(= 100) and z, normalised by theoretical 
prediction as 



D 



P(k;z)/P(k;z init ) 



Dp r cd -fpred(^) ^)/^~pred(&j ^init) 



(6) 



We compute the theoretical prediction directly from the re- 
sult of the Boltzmann solver. Run C and Run D reproduce 
the theoretically predicted value very well, whereas other 
runs show substantial deviations from the prediction. This 
remains true even at z = 30, as shown in Figure 3. 

The significant over-shoot of the baryonic component 
in Run A is explained by the tight coupling of the gas and 
CDM particles falsely caused by using correct transfer func- 
tions. In Figure 4, we plot the initial particle distributions in 
a slab of thickness 40 /i _1 kpc for Run A (top panel) and Run 
B (bottom panel). Figure 4 clearly shows that the gas par- 
ticles are placed too closely to the nearest CDM particles in 
Run A. Using separate transfer functions whose amplitudes 
differ by an approximately constant (and small) factor over 
the relevant length scales, one obtains a distribution of gas 
particles each of which is only slightly separated from the 
closest CDM particle. In this situation, a major fraction of 
the force exerted on a gas particle is due to the nearest CDM 
particle, which causes the gas particle to move faster than 
expected in linear theory. This false coupling was success- 
fully avoided in RunB, for which we used two sets of inde- 
pendent glass distributions; one for the gas particles and the 
other for the CDM particles. Thus, the velocity vector of a 
gas particle is uncorrelated with the direction to the nearest 
CDM particle, as is seen in Figure 4. 

Although Run B appears to behave reasonably well, 
CDM density fluctuations grow slightly too fast whereas the 
baryon density fluctuations grow too slowly, as seen in the 
growth rate in comparison with the theoretical prediction. 
The initial velocities are assigned separately to each com- 
ponent following the standard procedure of the Zel'dovich 
approximation. Without correcting velocities, the CDM den- 
sity fluctuations grow faster than the correct solution, and 
the opposite is true for the baryons. This feature is not seen 
in Run C, because we corrected the initial velocities for both 
the baryonic and CDM components. Indeed, the density fluc- 
tuations in Run C grow almost precisely as predicted from 
z init = 100 to z = 20. 

Run D reproduces the correct growth of density fluc- 
tuations. This is as expected, because, in Run D, the two 
components initially behave as a single component, until 
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Figure 4. The initial particle distribution in a slab of thickness 
40 h~ 1 kpc (1 percent of the box side-length) in Run A (top) and 
Run B (bottom). The solid points are CDM particles and the 
open diamonds are gas particles. The short lines originating from 
the gas particles show the velocity vectors. 



the pressure forces on gas particles becomes important and, 
making the orbits of the gas and CDM particles deviate 
from each other. Note that, at any output redshift, the den- 
sity fluctuations for a single component (baryons or CDM) 
does not accurately represent the true solution. Hence, we 
conclude that the set-up for Run D is not suitable for studies 
that are concerned with the difference in the distribution of 
baryons and dark matter. 

Overall the initial set up for Run C is clearly the most 
favourable. The method used for Run C appears to eliminate 
all the problems and undesirable features found in the other 
simulations. Also, the excellent agreement between the mea- 
sured power spectra and the theoretical prediction for both 
baryonic and CDM components over a large redshift interval 
confirms that the method is useful for the multi-component 
cosmological simulations. 
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Figure 5. The two-point correlation function for the gas par- 
ticles (£t>bi open diamonds) and for the dark matter particles 
(See, solid circles) in a glass distribution. The solid line shows the 
cross-correlation £bc computed for an arbitrary superposition of 
two independent glass files, whereas the dashed line is the cross- 
correlation measured for the distribution of the mixture of the two 
components after they were evolved for a brief period (see text). 
The separation length r is normalised by the mean inter-particle 
separation. 

5 DISCUSSION 

5.1 Glass distribution for two components 

Although the comparison of the results of Run A and Run B 
show the advantage of using two glass distributions for the 
CDM and baryonic components, it remains unclear whether 
the mixture of two independent glass distributions still pos- 
sesses all desired features. Generally it does not, because 
close encounters between particles from the different com- 
ponents could happen (Jun Makino, private communica- 
tion). In order to address this issue, we compute the cross- 
correlation function £bc for a mixed glass distribution. We 
make a distribution for 2 x 128 3 particles by arbitrarily 
superposing a glass distribution for 128 3 particles on an- 
other glass distribution which is independently generated. 
The solid line in Figure 5 shows the gas-dark matter cross- 
correlation £bc for this arbitrary mixture of two sets of glass 
distributions. The overall behaviour might seem quite good, 
as the cross-correlation fbc stays at an almost constant value 
around zero. However, it means that there is some chance 
that a gas particle finds a CDM particle quite close to it. 
Since this chance coupling of the two different components 
can cause exactly the same problem, even locally, as we 
found in Run A, it is preferable to generate a distribution of 
the two components such that the cross-correlation between 
them is reduced on very small scales. We follow the same 
procedure as we have done to create glass distributions; i.e., 
we switch the sign of the gravitational force and evolve the 
system consisting of a mixture of 2 x 128 3 particles in an ex- 
panding background. Since our primary purpose is to avoid 
false coupling of gas and dark matter particles, we evolve 
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Figure 6. The probablity distribution P(S) for the initial condi- 
tions for Run C (top) and Run D (bottom). In Run D, the CDM 
and gas particles have identical distibutions and hence P(5) is 
also identical initially. 



the system only for a brief period of time. Figure 5 shows 
the cross-correlation function for the initial (solid line) and 
evolved (dashed line) mixture of the particles. We also plot 
the correlation of the gas particles (open diamonds) and that 
of the CDM particles (solid circles) . Figure 5 shows that the 
glass distribution we obtained in this manner possesses the 
following desired properties: (1) the correlation of a single 
component is reduced on scales below the mean-interparticle 
separation, and (2) the cross-correlation is also reduced on 
the small scales, confirming that the false coupling with the 
other component such as that found in Run A will not be 
induced even when we use separate transfer functions. 
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Figure 7. The evolution of the two-point correlation function 
for Run C (top) and Run D (bottom). Symbols are for the dark 
matter component, as indicated in the legend in each panel, and 
the dashed line is the correlation function of the baryonic particles 
at 2 = 20. The separation length r is normalized by the mean 
inter-particle separation. 



5.2 Real space correlation on small scales 

While we conclude that glass distributions are preferred over 
a regular grid distribution for the reasons described in the 
previous sections, there still remain a few issues to be clar- 
ified regarding unperturbed particle distributions. Setting- 
up initial conditions for cosmological particle simulations is 
generally divided into two parts. One is the generation of a 
uniform particle distribution, which we have discussed, and 
the other is imposing desired density fluctuations by displac- 
ing the particles using the Zeldovich approximation (section 
2.2). For the currently standard cosmological models, we 
further assume that the initial density field is a Gaussian 



random field, for which the power spectrum fully describes 
the statistical properties. We first check and show that the 
initial particle distribution generated in this manner repre- 
sents a Gaussian density field, i.e., the probability distri- 
bution P(S) of overdensity <5(x) = p/p — 1 in real space is 
Gaussian. Figure 6 shows P(8) for our Run C and Run D. 
We compute the real space density 8 (x) by first re-sampling 
the particles onto a 128 3 grid and then smoothing on a scale 
R s = 100/i _1 kpc. For Run C, we compute P{8) for both 
the baryonic and dark matter components, whereas for Run 
D only the dark matter component is used. (Note that in 
Run D the gas particles are put on top of the dark mat- 
ter particles initially, and hence P(8) at the initial epoch is 
identical.) As clearly seen in Figure 6, the computed proba- 
bility distribution is well fitted by a Gaussian for each case. 
The difference between the two components in Run C is due 
to the difference in the linear power amplitude on the cho- 
sen length scale R s ~ 100/i~ 1 kpc. The measured variances 
for the baryonic and dark matter components are <7b=0.062, 
C7cdm=0.081, respectively, in good agreement with analytical 
estimates a thib =0.060, o" thjCdm =0.083 computed from the in- 
put density power spectra. The corresponding numbers for 
Run D are a to tai=0.076 and cr t h,totai=0.078. We note that, 
in practice, we used a random number generator in assign- 
ing the phases for each Fourier mode so that the generated 
initial conditions possess a desired feature of a Gaussian 
random field. 

Since any particle distribution (e.g. glass or grid) has 
its own intrinsic correlations (Gabrielli et al. 2002), the 
perturbed distribution will have, in principle, a superpo- 
sition of the intrinsic and imposed (desired) correlations. 
Baertschiger & Sylos-Labini (2002) argue that the desired 
power law correlations for cosmological models are not prop- 
erly realized in discritised density fields at the initial epoch 
(see also Knebe & Dominguez 2003) . We quantify and show 
this feature by computing the real space correlation from 
particle distributions. We use the pair-count estimator pro- 
posed by Landy & Szalay (1993): 



£(r) = 



DD(r) - 2DR(r) + RR(r) 
RR(r) 



(7) 



We randomly distribute the same number of particles as the 
actual simulation particle, and evaluate the data-data (DD), 
data-random (DR), and random-random (RR) terms to ob- 
tain £(r). In Figure 7 we show the two-point correlation func- 
tions measured from the outputs of Run C and Run D. Ini- 
tially the correlation on scales below the mean inter-particle 
separation (dotted line) is absent, reflecting the damped cor- 
relation on the small scales in the unperturbed distribution 
(see Figure 5) . The particle clustering at early epochs on the 
smallest scales is considerably affected by the discreteness as 
argued by Hamana, Yoshida & Suto (2002). Baertschiger et 
al. (2002) claim that the development of a power-law corre- 
lation on the small scales is governed partly by the initial 
clustering. It is interesting that the two-point correlation 
functions of the two components at z = 20 are quite simi- 
lar even on the smallest length scales, despite the fact that 
they started from different glass distributions. On the other 
hand, the feature that the correlation amplitudes are differ- 
ent between the two components in Run C is plausible, and 
it appears to reproduce correctly the smaller density fluctu- 
ations of baryons. At late epochs the nonlinear power trans- 
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fer eventually dominates the initial power if the effective 
spectral index at the scale n = dlog P(fc)/dlog k is much 
less than -1. CDM models satisfy this condition on small 
length scales (Suto & Sasaki f99f), and hence the power on 
scales below the mean particle separation at later epochs is 
expected to be generated via nonlinear mode-coupling (Sug- 
inohara et al. 199f ; Hamana et al. 2002). Further theoretical 
study is clearly needed to understand the behaviour of the 
small scale correlations at early epochs. 

5.3 Differing phases 

As we discussed in Section 3.3, the phases of the pertur- 
bations are assumed to be identical for baryons and dark 
matter as long as the density perturbations are in the lin- 
ear regime. Nevertheless, it is interesting to examine how 
the growth of density fluctuations is affected by offsets of 
phases between the two components. We note that, in real 
astrophysical situations, phase offsets may be caused during 
early reionization when the intergalactic gas is photo-heated 
by radiation from an early generation of stars and driven out 
of dark matter potential wells. In such cases, it is expected 
that the growth of either or both components is delayed, 
compared to the case where the phases are synchronised. 
We carry out a simple test by applying independent phases 
to baryon and CDM perturbations at the initial epoch. For 
simplicity, we assign completely random phases, rather than 
specifying actual mechanisms that could cause such an off- 
set. Hence, the test is meant to show the effects of the phase 
offset qualitatively. Figure 8 shows the evolution of the den- 
sity fluctuations for Run E. As expected, the evolution of 
both the baryon and the CDM density fluctuations is de- 
layed, with the effect on the former being more substan- 
tial. On the other hand, it is interesting that the growth 
of the baryon density fluctuations is not entirely prevented, 
despite the large difference in the fractional densities be- 
tween the two components. This is because, on large scales, 
the dark matter density fluctuations are still small and the 
gravitational force from the dark matter has not significantly 
changed the initial momentum of the baryonic components. 
Note that, in Figure 8, the apparent differences between 
the measured power spectra and the theoretical prediction 
should be interpreted to be owing mostly to the difference in 
the background assumption, because the theoretical predic- 
tion itself is obtained assuming the initial phases are identi- 
cal. 

5.4 Nonlinear objects 

In the previous sections, we have focused on the evolution of 
density fluctuations on large scales when they are in linear 
regime. As Figure 2 and Figure 3 show, the overall differ- 
ences between our simulations appear as if their effective 
growth factors deviate, falsely or due to an inaccurate set- 
up, from the true value. This raises an interesting question 
of whether or not such differences also show up in some 
properties of non-linear objects. Using our simulations, we 
identify dark matter halos by running a friend of friends 
(FOF) group-finder with the standard choice of linking pa- 
rameter b — 0.164. We first consider only the dark matter 
component. The baryon fraction within the halos will be ex- 
amined later in this section. Figure 9 shows the cumulative 
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Figure 8. The evolution of the power spectra in Run E. As in 
Figure 3, the power spectra at z = 20 are scaled vertically by a 
constant factor of five, in order to clarify the plot. 



mass functions for our simulations at z — 20 and z = 15. 
Since the overall abundance of halos is still small at z = 20, 
as seen in the top panel in Figure 9, we allow our simula- 
tions to evolve to z — 15 to have a large number of halos. 
We show the result for the z = 15 outputs in the bottom 
panel of Figure 9. We compute the mass of the identified ha- 
los to be the number of member particles times the particle 
mass. We discard halos consisting of fewer than 25 parti- 
cles, and thus the minimum mass in our sample is set to be 
5 x 10 7 hT 1 Mq. Although the mass functions measured in 
the four simulations agree reasonably well, there is an ap- 
preciable differences, particularly between Run D and the 
other runs. This is easily accounted for by the fact that the 
growth of the dark matter density fluctuations in Run D is 
slower than the other runs (see Figure 2, 3). 

We also measured the baryon fraction within the dark 
halos. To this end, we compute the virial radii and virial 
masses for the halos following Yoshida et al. (2002) . Briefly, 
we define the virial radius R v i T as the radius of the sphere 
centred on the most bound particle of the FOF group hav- 
ing overdensity 200 with respect to the critical density. The 
virial mass Af v i r is then the enclosed mass (gas and dark 
matter) within R v i r . We compute the baryon fractions for 
the halos as /b = M gas /Af total, and normalise by the global 
baryon fraction fib /fimatter- Figure 10 shows the measured 
baryon fractions for the halos in our simulation at z = 15. 
In all cases the (normalised) baryon fractions are smaller 
than unity, representing correctly the effect of hydrodynamic 
pressure in the absence of gas cooling. Whereas the large 
scatter at masses smaller than ~ 2 x 10 s /i _1 M0, consisting 
of fewer than 100 particles, is due to the mass resolution 
(see Yoshida et al. 2002 for a discussion of the resolution 
issue), the difference on larger mass scale between the four 
runs is unexpected. Again, by noting the residual difference 
in the baryon and dark matter density fluctuations shown 
in Figure 3 and Figure 7, we are led to the conclusion that 
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Figure 9. Cumulative mass function of dark matter halos at 
2 = 20 (top) and z = 15 (bottom). The vertical lines indicate the 
number of particles constituting a halo of a corresponding mass. 



the difference in the density fluctuations at early epochs is 
responsible for the overall trend shown in Figure 10. Indeed, 
the amplitudes of the baryon and dark matter density fluc- 
tuations in Run A and Run D are quite similar from early 
on (z ~ 50), while those in Run B and Run C remain rela- 
tively large until low redshift (z ~ 20). Thus the gas infall 
is slightly delayed in Run B and Run C, resulting in the 
smaller baryon fractions than in Run A and Run D. 

It should be noted that these results are not expected 
to be precise in details, but are likely to be dependent on the 
procedure for identifying halos and also on the definition of 
the halo mass. The detailed description of nonlinear evolu- 
tion is beyond the scope of our paper and thus we reserve this 
issue for the subject of future work using simulations with 
substantially higher resolution than those described here. 
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Figure 10. Baryon fractions for massive (M > 10 8 /i -1 Mq) halos 
at 2 = 15. The baryon fractions are normalised by the global 
baryon fraction fi baryon /S! ma ttcr- 

6 CONCLUSION 

We have studied the evolution of the matter density fluctu- 
ations in cosmological A^-body/hydrodynamic simulations, 
and discussed various aspects of the set-up of the initial 
conditions. All the methods we used in generating the ini- 
tial conditions except for Run E may appear to be plausible, 
differing only in a few details. However, we have found that 
such seemingly minor details cause appreciable differences 
in the evolved density field, not only in linear regime but 
also in the non-linear regime. 

Interestingly, it has turned out that neither of the meth- 
ods used for Run A and Run D provide suitable initial con- 
ditions, with the results of the former being more signifi- 
cantly affected by numerical artifacts. This perhaps comes 
as a surprise, particularly in the case of Run A, because one 
may naively hope that using separate transfer functions for 
baryons and CDM would produce more accurate results. In 
fact, the baryon density fluctuations reproduced in Run A 
are the least accurate in terms of the growth rate. Although 
the particular problem of false particle coupling found in 
Run A is likely to be a problem in simulations where both 
components are represented by discrete mass elements, it 
is preferable, even in simulations that employ a grid-based 
Eulerian scheme, to avoid artificial coherence between the 
initial fluid velocity vector and the direction to the near- 
est CDM particle. This may deserve further study using a 
grid-based scheme. 

By analysing the outputs of the numerical experi- 
ments, we identified some peculiar problems and inaccura- 
cies caused by the simplest methods. Then, based on the re- 
sults, we have suggested a more accurate method to set up 
initial conditions; i.e., the one used for Run C. In summary, 
we recommend: (1) using independent glass particle distri- 
butions from which a "glass mixture" is made by a treatment 
to avoid false coupling, (2) using two transfer functions com- 
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puted for baryons and for dark matter, and (3) taking into 
account the difference in their velocities at the initialisation 
epoch. We emphasize that the latter two procedures may be 
necessary for all cosmological simulations employing multi- 
ple components, whether or not both dark matter and gas 
are realised with particles. We have explicitly shown that, 
in this manner, the evolution of matter density fluctuations 
are accurately followed at high redshift. 
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